clear all
set more off

global data 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_dta"
global figures 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_log"
global in 		"R:\SharedProjects\Shared2020-070\2016\input"

cap log close
log using $figures\Q_file,replace t

***Run this first for order_expansion=3 and then for order_expansion=2 (line 1048)*/
********************************************************************************
scalar define order_expansion=3		/*2=quadratic expansion, 3=cubic expansion*/
********************************************************************************

cd $data

************************************************************************************************************************************************************************************
******VIGNETTE DATA
************************************************************************************************************************************************************************************
************************************************************************************************************************************************************************************
u $data\vignette_data, clear
******************************************************************************************************************************************************************

ren walkra walkra_R_1 
ren dressa dressa_R_1 
ren stoopa stoopa_R_1 
ren beda   beda_R_1   
ren bmi    bmi_R      
ren hosp   hosp_R      

qui tab occ_interm,gen(occx)
qui tab vign_id,gen(vigndomx)

************************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global spec_L 		female college racex2 doctor* married widowed age  walkra dressa stoopa beda hosp bmi i.occ_interm experience
global spec_V 		$spec_L vign_f i.vign_id 
************************************************************************************************************************************************************************************

******************************************************************************************************************************************************************
**Sample selection (drop those w/ missing data in the variables used for the 2 regressions)
qui biprobit (vign_severe = $spec_V) (cantwork = $spec_L) if year==2007, vce(cluster hhidpn)	
keep if e(sample)
******************************************************************************************************************************************************************

*This is the bivariate probit L(respondent), L(vignette)
****TABLE 9 + \RHO_{L,V}
eststo bigml2: biprobit (vign_severe = $spec_V) (cantwork = $spec_L) if year==2007, vce(cluster hhidpn)	

matrix bV=e(b)
matrix bV=bV'

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments	=_b[vign_severe:female]\_b[vign_severe:vign_f]\r(table)[1,colsof(r(table))]
matrix semoments=_se[vign_severe:female]\_se[vign_severe:vign_f]\r(table)[2,colsof(r(table))]	
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
************************************************************************************************************************************************

************************************************************************************************************************************************
estadd scalar corr_L_V 	r(table)[1,colsof(r(table))]


************************************************************************************************************************************************************************************
********REJECTION DATA *************************************************************************************************************************************************************************************
*************************************************************************************************************************************************************************************

u $data\rejection_data, clear
***QUADRATIC EXPANSION
g agesq=age^2
g expsq=experience^2

foreach x of varlist age experience agesq expsq college racex2 married widowed walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R {
	gen `x'fem=`x'*female
}

*****************************
if order_expansion==3	{	
	gen agecb=age^3
	gen expcb=experience^3
	gen agecbfem=agecb*female
	gen expcbfem=expcb*female
}
*****************************

forvalues j=1(1)9 {
	gen bsx`j'fem=bsx`j'*female
}
forvalues j=1(1)8 {
	gen occint`j'fem=(occ_interm==`j')*female
}
save $data\rejection_data_extended,replace

u $data\rejection_data_extended, clear

/*This is -X*a(Lbar)*/
#delimit ;
gen a1_fromV=
bV[2 ,1]		*college+  
bV[3 ,1]         *racex2  +
bV[4 ,1]   *doctor_told_has_hbp +
bV[5 ,1]   *doctor_told_has_psy +
bV[6 ,1]   *doctor_told_has_hea +
bV[7 ,1]   *doctor_told_has_art +
bV[8 ,1]   *doctor_told_has_dia +
bV[9 ,1]   *doctor_told_has_lun +
bV[10,1]   *doctor_told_has_str +
bV[11,1]   *doctor_told_has_can +
bV[12,1]         *married +
bV[13,1]         *widowed +
bV[14,1]             *age +
bV[15,1]      *walkra_R_1 +
bV[16,1]      *dressa_R_1 +
bV[17,1]      *stoopa_R_1 +
bV[18,1]        *beda_R_1 +
bV[19,1]          *hosp_R +
bV[20,1]          *bmi_R  +
bV[21,1]     *(occ_interm==1)+
bV[22,1]     *(occ_interm==2)+
bV[23,1]     *(occ_interm==3)+
bV[24,1]     *(occ_interm==4)+
bV[25,1]     *(occ_interm==5)+
bV[26,1]     *(occ_interm==6)+
bV[27,1]     *(occ_interm==7)+
bV[28,1]     *(occ_interm==8)+
bV[29,1]      *experience;
#delimit cr

/*This is -X*a(Lbar)(institutional variables only)*/
#delimit ;
gen a1_inst_fromV=
bV[2 ,1]		*college+  
bV[4 ,1]   *doctor_told_has_hbp +
bV[5 ,1]   *doctor_told_has_psy +
bV[6 ,1]   *doctor_told_has_hea +
bV[7 ,1]   *doctor_told_has_art +
bV[8 ,1]   *doctor_told_has_dia +
bV[9 ,1]   *doctor_told_has_lun +
bV[10,1]   *doctor_told_has_str +
bV[11,1]   *doctor_told_has_can +
bV[14,1]             *age +
bV[15,1]      *walkra_R_1 +
bV[16,1]      *dressa_R_1 +
bV[17,1]      *stoopa_R_1 +
bV[18,1]        *beda_R_1 +
bV[19,1]          *hosp_R +
bV[20,1]          *bmi_R  +
bV[21,1]     *(occ_interm==1)+
bV[22,1]     *(occ_interm==2)+
bV[23,1]     *(occ_interm==3)+
bV[24,1]     *(occ_interm==4)+
bV[25,1]     *(occ_interm==5)+
bV[26,1]     *(occ_interm==6)+
bV[27,1]     *(occ_interm==7)+
bV[28,1]     *(occ_interm==8)+
bV[29,1]      *experience;
#delimit cr


************************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global c_L 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm 
global c_A 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm logbs 
global c_R 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx1-bsx9 agesq-occint8fem 
************************************************************************************************************************************************************************************

gen A=applied
gen R=nosuccess
gen L=cantwork_notemp

gen 	index2=1 if R==1 & A==1
replace index2=2 if R==0 & A==1
replace index2=3 if A==0

set seed 2134569

// Create Halton draw variables
mdraws, dr(10) neq(3) prefix(z)			/*If dr(N), then in egen `sp3' and egen `sp2' adjust dr(N) accordingly*/

// Get initial values
quietly {
        probit A $c_A
        mat b1 = e(b)
        mat coleq b1 = A
        probit L $c_L
        mat b2 = e(b)
        mat coleq b2 = L
        probit R $c_R
        mat b3 = e(b)
        mat coleq b3 = R
        mat b0 = b1, b2, b3
}


cap program drop myll2
program define myll2
args lnf xb1 xb2 xb3 c21 c31 c32
        tempvar sp2 sp3 k1 k2 k3
        quietly {
            gen double `k1' = 2*$ML_y1 - 1
            gen double `k2' = 2*$ML_y2 - 1
            gen double `k3' = 2*$ML_y3 - 1
            tempname cf21 cf31 cf32 cf22 cf33 C1 C2
            su `c21', meanonly
            scalar `cf21' = r(mean)
            su `c31', meanonly
            scalar `cf31' = r(mean)
            su `c32', meanonly
            scalar `cf32' = r(mean)
            scalar `cf22' = sqrt( 1 - `cf21'^2 )
            scalar `cf33' = sqrt( 1 - `cf31'^2 - `cf32'^2 )
            tempname C1 C2
            mat `C1' = (1, 0, 0 \ `cf21', `cf22', 0 \ `cf31', `cf32', `cf33')
            mat `C2' = (1, 0 \ `cf21', `cf22' )
            egen `sp3' = mvnp(`xb1' `xb2' `xb3') if $ML_y1==1, ///
                    chol(`C1') dr(10) prefix(z) signs(`k1' `k2' `k3') adoonly
            egen `sp2' = mvnp(`xb1' `xb2' ) if $ML_y1==0, ///
                    chol(`C2') dr(10) prefix(z) signs(`k1' `k2') adoonly
            replace `lnf'= ln(`sp3') if $ML_y1==1
            replace `lnf'= ln(`sp2') if $ML_y1==0
        }
end

ml model lf myll2 	(A: A=$c_A) ///
					(L: L=$c_L) ///
					(R: R=$c_R) ///
					/c21 /c31 /c32 ///
					, missing

ml init b0
eststo bigml1:ml maximize

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[L:female]\_b[A:female]
matrix semoments=semoments\_se[L:female]\_se[A:female]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar betaFL=_b[L:female]	
scalar delta =_b[R:female]	/*to be used for pinning down \alpha*/

gen bF_L=_b[L:female]
gen bF_A=_b[A:female] 
gen bF_R=_b[R:female]
	
#delimit;
gen a1_age=_b[L:age]*age;
gen a1_institutional_year=
_b[L:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[L:doctor_told_has_psy]* doctor_told_has_psy+
_b[L:doctor_told_has_hea]* doctor_told_has_hea+
_b[L:doctor_told_has_art]* doctor_told_has_art+
_b[L:doctor_told_has_dia]* doctor_told_has_dia+
_b[L:doctor_told_has_lun]* doctor_told_has_lun+
_b[L:doctor_told_has_str]* doctor_told_has_str+
_b[L:doctor_told_has_can]* doctor_told_has_can+
_b[L:walkra_R_1]*             walkra_R_1+
_b[L:dressa_R_1]*             dressa_R_1+
_b[L:stoopa_R_1]*             stoopa_R_1+
_b[L:beda_R_1]*                 beda_R_1+
_b[L:bmi_R]*                       bmi_R+
_b[L:hosp_R]*                     hosp_R+
_b[L:age]*                           age+
_b[L:college]*                   college+
_b[L:experience]*             experience+
_b[L:1b.occ_interm]*        (occ_interm==1)+
_b[L:2.occ_interm]*         (occ_interm==2)+
_b[L:3.occ_interm]*         (occ_interm==3)+
_b[L:4.occ_interm]*         (occ_interm==4)+
_b[L:5.occ_interm]*         (occ_interm==5)+
_b[L:6.occ_interm]*         (occ_interm==6)+
_b[L:7.occ_interm]*         (occ_interm==7)+
_b[L:8.occ_interm]*         (occ_interm==8)+
_b[L:1b.wave]*                  (wave==1)+
_b[L:2.wave]*                   (wave==2)+
_b[L:3.wave]*                   (wave==3)+
_b[L:4.wave]*                   (wave==4)+
_b[L:5.wave]*                   (wave==5)+
_b[L:6.wave]*                   (wave==6)+
_b[L:7.wave]*                   (wave==7)+
_b[L:8.wave]*                   (wave==8)+
_b[L:9.wave]*                   (wave==9)+
_b[L:10.wave]*                  (wave==10)+
_b[L:11.wave]*                  (wave==11)+
_b[L:12.wave]*                  (wave==12)+
_b[L:13.wave]*                  (wave==13)+
_b[L:14.wave]*                  (wave==14)+
_b[L:15.wave]*                  (wave==15);
gen a1_health=
_b[L:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[L:doctor_told_has_psy]* doctor_told_has_psy+
_b[L:doctor_told_has_hea]* doctor_told_has_hea+
_b[L:doctor_told_has_art]* doctor_told_has_art+
_b[L:doctor_told_has_dia]* doctor_told_has_dia+
_b[L:doctor_told_has_lun]* doctor_told_has_lun+
_b[L:doctor_told_has_str]* doctor_told_has_str+
_b[L:doctor_told_has_can]* doctor_told_has_can+
_b[L:walkra_R_1]*             walkra_R_1+
_b[L:dressa_R_1]*             dressa_R_1+
_b[L:stoopa_R_1]*             stoopa_R_1+
_b[L:beda_R_1]*                 beda_R_1+
_b[L:bmi_R]*                       bmi_R+
_b[L:hosp_R]*                     hosp_R;
gen a1_demo=
_b[L:female]*				      female+
_b[L:age]*                           age+
_b[L:racex2]*                     racex2+
_b[L:married]*                   married+
_b[L:widowed]*                   widowed+
_b[L:_cons];
gen a1_labor=
_b[L:college]*                   college+
_b[L:experience]*             experience+
_b[L:1b.wave]*                  (wave==1)+
_b[L:2.wave]*                   (wave==2)+
_b[L:3.wave]*                   (wave==3)+
_b[L:4.wave]*                   (wave==4)+
_b[L:5.wave]*                   (wave==5)+
_b[L:6.wave]*                   (wave==6)+
_b[L:7.wave]*                   (wave==7)+
_b[L:8.wave]*                   (wave==8)+
_b[L:9.wave]*                   (wave==9)+
_b[L:10.wave]*                  (wave==10)+
_b[L:11.wave]*                  (wave==11)+
_b[L:12.wave]*                  (wave==12)+
_b[L:13.wave]*                  (wave==13)+
_b[L:14.wave]*                  (wave==14)+
_b[L:15.wave]*                  (wave==15)+
_b[L:1b.occ_interm]*        (occ_interm==1)+
_b[L:2.occ_interm]*         (occ_interm==2)+
_b[L:3.occ_interm]*         (occ_interm==3)+
_b[L:4.occ_interm]*         (occ_interm==4)+
_b[L:5.occ_interm]*         (occ_interm==5)+
_b[L:6.occ_interm]*         (occ_interm==6)+
_b[L:7.occ_interm]*         (occ_interm==7)+
_b[L:8.occ_interm]*         (occ_interm==8);
#delimit cr
gen 	a1=a1_demo+a1_health+a1_labor
foreach x in a1 a1_institutional_year a1_demo a1_health a1_labor a1_age	{
	replace `x'=. if L==.
}

#delimit;
gen a2_age=_b[A:age]*age;
gen a2_institutional_year=
_b[A:age]*                           age+
_b[A:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[A:doctor_told_has_psy]* doctor_told_has_psy+
_b[A:doctor_told_has_hea]* doctor_told_has_hea+
_b[A:doctor_told_has_art]* doctor_told_has_art+
_b[A:doctor_told_has_dia]* doctor_told_has_dia+
_b[A:doctor_told_has_lun]* doctor_told_has_lun+
_b[A:doctor_told_has_str]* doctor_told_has_str+
_b[A:doctor_told_has_can]* doctor_told_has_can+
_b[A:walkra_R_1]*             walkra_R_1+
_b[A:dressa_R_1]*             dressa_R_1+
_b[A:stoopa_R_1]*             stoopa_R_1+
_b[A:beda_R_1]*                 beda_R_1+
_b[A:bmi_R]*                       bmi_R+
_b[A:hosp_R]*                     hosp_R+
_b[A:college]*                   college+
_b[A:experience]*             experience+
_b[A:1b.occ_interm]*        (occ_interm==1)+
_b[A:2.occ_interm]*         (occ_interm==2)+
_b[A:3.occ_interm]*         (occ_interm==3)+
_b[A:4.occ_interm]*         (occ_interm==4)+
_b[A:5.occ_interm]*         (occ_interm==5)+
_b[A:6.occ_interm]*         (occ_interm==6)+
_b[A:7.occ_interm]*         (occ_interm==7)+
_b[A:8.occ_interm]*         (occ_interm==8)+
_b[A:1b.wave]*                  (wave==1)+
_b[A:2.wave]*                   (wave==2)+
_b[A:3.wave]*                   (wave==3)+
_b[A:4.wave]*                   (wave==4)+
_b[A:5.wave]*                   (wave==5)+
_b[A:6.wave]*                   (wave==6)+
_b[A:7.wave]*                   (wave==7)+
_b[A:8.wave]*                   (wave==8)+
_b[A:9.wave]*                   (wave==9)+
_b[A:10.wave]*                  (wave==10)+
_b[A:11.wave]*                  (wave==11)+
_b[A:12.wave]*                  (wave==12)+
_b[A:13.wave]*                  (wave==13)+
_b[A:14.wave]*                  (wave==14)+
_b[A:15.wave]*                  (wave==15);
gen a2_demo=
_b[A:female]*				      female+
_b[A:age]*                           age+
_b[A:racex2]*                     racex2+
_b[A:married]*                   married+
_b[A:widowed]*                   widowed+
_b[A:logbs]*logbs						+
_b[A:_cons];
gen a2_health=
_b[A:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[A:doctor_told_has_psy]* doctor_told_has_psy+
_b[A:doctor_told_has_hea]* doctor_told_has_hea+
_b[A:doctor_told_has_art]* doctor_told_has_art+
_b[A:doctor_told_has_dia]* doctor_told_has_dia+
_b[A:doctor_told_has_lun]* doctor_told_has_lun+
_b[A:doctor_told_has_str]* doctor_told_has_str+
_b[A:doctor_told_has_can]* doctor_told_has_can+
_b[A:walkra_R_1]*             walkra_R_1+
_b[A:dressa_R_1]*             dressa_R_1+
_b[A:stoopa_R_1]*             stoopa_R_1+
_b[A:beda_R_1]*                 beda_R_1+
_b[A:bmi_R]*                       bmi_R+
_b[A:hosp_R]*                     hosp_R;
gen a2_labor=
_b[A:college]*                   college+
_b[A:experience]*             experience+
_b[A:1b.wave]*                  (wave==1)+
_b[A:2.wave]*                   (wave==2)+
_b[A:3.wave]*                   (wave==3)+
_b[A:4.wave]*                   (wave==4)+
_b[A:5.wave]*                   (wave==5)+
_b[A:6.wave]*                   (wave==6)+
_b[A:7.wave]*                   (wave==7)+
_b[A:8.wave]*                   (wave==8)+
_b[A:9.wave]*                   (wave==9)+
_b[A:10.wave]*                  (wave==10)+
_b[A:11.wave]*                  (wave==11)+
_b[A:12.wave]*                  (wave==12)+
_b[A:13.wave]*                  (wave==13)+
_b[A:14.wave]*                  (wave==14)+
_b[A:15.wave]*                  (wave==15)+
_b[A:1b.occ_interm]*        (occ_interm==1)+
_b[A:2.occ_interm]*         (occ_interm==2)+
_b[A:3.occ_interm]*         (occ_interm==3)+
_b[A:4.occ_interm]*         (occ_interm==4)+
_b[A:5.occ_interm]*         (occ_interm==5)+
_b[A:6.occ_interm]*         (occ_interm==6)+
_b[A:7.occ_interm]*         (occ_interm==7)+
_b[A:8.occ_interm]*         (occ_interm==8);
#delimit cr
gen 	a2=a2_demo+a2_health+a2_labor
foreach x in a2 a2_institutional_year a2_demo a2_health a2_labor a2_age	{
	replace `x'=. if A==.
}

#delimit;
gen a3_age=_b[R:age]*age;
gen a3_institutional_year=
	_b[R:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[R:doctor_told_has_psy]* doctor_told_has_psy+
	_b[R:doctor_told_has_hea]* doctor_told_has_hea+
	_b[R:doctor_told_has_art]* doctor_told_has_art+
	_b[R:doctor_told_has_dia]* doctor_told_has_dia+
	_b[R:doctor_told_has_lun]* doctor_told_has_lun+
	_b[R:doctor_told_has_str]* doctor_told_has_str+
	_b[R:doctor_told_has_can]* doctor_told_has_can+
	_b[R:walkra_R_1]*             walkra_R_1+
	_b[R:dressa_R_1]*             dressa_R_1+
	_b[R:stoopa_R_1]*             stoopa_R_1+
	_b[R:beda_R_1]*                 beda_R_1+
	_b[R:bmi_R]*                       bmi_R+
	_b[R:hosp_R]*                     hosp_R+
	_b[R:bsx1]*bsx1+
	_b[R:bsx2]*bsx2+
	_b[R:bsx3]*bsx3+
	_b[R:bsx4]*bsx4+
	_b[R:bsx5]*bsx5+
	_b[R:bsx6]*bsx6+
	_b[R:bsx7]*bsx7+
	_b[R:bsx8]*bsx8+
	_b[R:bsx9]*bsx9+
	_b[R:age]*                           age+
	_b[R:college]*                   college+
	_b[R:experience]*             experience+
	_b[R:1b.occ_interm]*        (occ_interm==1)+ 
	_b[R:2.occ_interm]*         (occ_interm==2)+
	_b[R:3.occ_interm]*         (occ_interm==3)+
	_b[R:4.occ_interm]*         (occ_interm==4)+
	_b[R:5.occ_interm]*         (occ_interm==5)+
	_b[R:6.occ_interm]*         (occ_interm==6)+
	_b[R:7.occ_interm]*         (occ_interm==7)+
	_b[R:8.occ_interm]*         (occ_interm==8)+
	_b[R:1991b.year]*(year==1991)+
	_b[R:1992.year] *(year==1992)+
	_b[R:1993.year] *(year==1993)+
	_b[R:1994.year] *(year==1994)+
	_b[R:1995.year] *(year==1995)+
	_b[R:1996.year] *(year==1996)+
	_b[R:1997.year] *(year==1997)+
	_b[R:1998.year] *(year==1998)+
	_b[R:1999.year] *(year==1999)+
	_b[R:2000.year] *(year==2000)+
	_b[R:2001.year] *(year==2001)+
	_b[R:2002.year] *(year==2002)+
	_b[R:2003.year] *(year==2003)+
	_b[R:2004.year] *(year==2004)+
	_b[R:2005.year] *(year==2005)+
	_b[R:2006.year] *(year==2006)+
	_b[R:2007.year] *(year==2007)+
	_b[R:2008.year] *(year==2008)+
	_b[R:2009.year] *(year==2009)+
	_b[R:2010.year] *(year==2010)+
	_b[R:2011.year] *(year==2011)+
	_b[R:2012.year] *(year==2012)+
	_b[R:2013.year] *(year==2013)+
	_b[R:2014.year] *(year==2014)+
	_b[R:2015.year] *(year==2015)+
	_b[R:2016.year] *(year==2016)+
	_b[R:2017.year] *(year==2017)+
	_b[R:2018.year] *(year==2018)+
	_b[R:2019.year] *(year==2019)+
	_b[R:2020.year] *(year==2020);	
gen a3_health=
	_b[R:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[R:doctor_told_has_psy]* doctor_told_has_psy+
	_b[R:doctor_told_has_hea]* doctor_told_has_hea+
	_b[R:doctor_told_has_art]* doctor_told_has_art+
	_b[R:doctor_told_has_dia]* doctor_told_has_dia+
	_b[R:doctor_told_has_lun]* doctor_told_has_lun+
	_b[R:doctor_told_has_str]* doctor_told_has_str+
	_b[R:doctor_told_has_can]* doctor_told_has_can+
	_b[R:walkra_R_1]*             walkra_R_1+
	_b[R:dressa_R_1]*             dressa_R_1+
	_b[R:stoopa_R_1]*             stoopa_R_1+
	_b[R:beda_R_1]*                 beda_R_1+
	_b[R:bmi_R]*                       bmi_R+
	_b[R:hosp_R]*                     hosp_R+
	_b[R:bsx1]*bsx1+
	_b[R:bsx2]*bsx2+
	_b[R:bsx3]*bsx3+
	_b[R:bsx4]*bsx4+
	_b[R:bsx5]*bsx5+
	_b[R:bsx6]*bsx6+
	_b[R:bsx7]*bsx7+
	_b[R:bsx8]*bsx8+
	_b[R:bsx9]*bsx9;
gen a3_demo=
	_b[R:female]*				      female+
	_b[R:age]*                           age+
	_b[R:racex2]*                     racex2+
	_b[R:married]*                   married+
	_b[R:widowed]*                   widowed+
	_b[R:_cons];
gen a3_labor=
	_b[R:college]*                   college+
	_b[R:experience]*             experience+
	_b[R:1991b.year]*(year==1991)+
	_b[R:1992.year] *(year==1992)+
	_b[R:1993.year] *(year==1993)+
	_b[R:1994.year] *(year==1994)+
	_b[R:1995.year] *(year==1995)+
	_b[R:1996.year] *(year==1996)+
	_b[R:1997.year] *(year==1997)+
	_b[R:1998.year] *(year==1998)+
	_b[R:1999.year] *(year==1999)+
	_b[R:2000.year] *(year==2000)+
	_b[R:2001.year] *(year==2001)+
	_b[R:2002.year] *(year==2002)+
	_b[R:2003.year] *(year==2003)+
	_b[R:2004.year] *(year==2004)+
	_b[R:2005.year] *(year==2005)+
	_b[R:2006.year] *(year==2006)+
	_b[R:2007.year] *(year==2007)+
	_b[R:2008.year] *(year==2008)+
	_b[R:2009.year] *(year==2009)+
	_b[R:2010.year] *(year==2010)+
	_b[R:2011.year] *(year==2011)+
	_b[R:2012.year] *(year==2012)+
	_b[R:2013.year] *(year==2013)+
	_b[R:2014.year] *(year==2014)+
	_b[R:2015.year] *(year==2015)+
	_b[R:2016.year] *(year==2016)+
	_b[R:2017.year] *(year==2017)+
	_b[R:2018.year] *(year==2018)+
	_b[R:2019.year] *(year==2019)+
	_b[R:2020.year] *(year==2020)+
	_b[R:1b.occ_interm]*        (occ_interm==1)+ 
	_b[R:2.occ_interm]*         (occ_interm==2)+
	_b[R:3.occ_interm]*         (occ_interm==3)+
	_b[R:4.occ_interm]*         (occ_interm==4)+
	_b[R:5.occ_interm]*         (occ_interm==5)+
	_b[R:6.occ_interm]*         (occ_interm==6)+
	_b[R:7.occ_interm]*         (occ_interm==7)+
	_b[R:8.occ_interm]*         (occ_interm==8);
#delimit cr

if order_expansion==2	{	
	gen a3_interaction= 					 	///
	_b[R:agesq] *          agesq			+	///	
	_b[R:expsq] *          expsq		    +	///
	_b[R:agefem] *     agefem		+			///
	_b[R:agesqfem] *     agesqfem		+		///
	_b[R:experiencefem] *     experiencefem	+   ///
	_b[R:expsqfem] *     expsqfem		+       ///
	_b[R:collegefem] *     collegefem		+   ///
	_b[R:racex2fem] *      racex2fem		+   ///
	_b[R:marriedfem] *     marriedfem		+   ///
	_b[R:widowedfem] *     widowedfem		+	///
	_b[R:walkra_R_1fem] *  walkra_R_1fem    +	///
	_b[R:dressa_R_1fem] *  dressa_R_1fem    +	///
	_b[R:stoopa_R_1fem] *  stoopa_R_1fem    +	///
	_b[R:beda_R_1fem] *    beda_R_1fem      +   ///
	_b[R:bmi_Rfem] *       bmi_Rfem         +   ///
	_b[R:hosp_Rfem] *      hosp_Rfem        +   ///
	_b[R:bsx1fem] *        bsx1fem          +   ///
	_b[R:bsx2fem] *        bsx2fem          +   ///
	_b[R:bsx3fem] *        bsx3fem          +   ///
	_b[R:bsx4fem] *        bsx4fem          +   ///
	_b[R:bsx5fem] *        bsx5fem          +   ///
	_b[R:bsx6fem] *        bsx6fem          +   ///
	_b[R:bsx7fem] *        bsx7fem          +   ///
	_b[R:bsx8fem] *        bsx8fem          +   ///
	_b[R:bsx9fem] *        bsx9fem          +   ///
	_b[R:occint1fem] *     occint1fem       +   ///
	_b[R:occint2fem] *     occint2fem       +   ///
	_b[R:occint3fem] *     occint3fem       +   ///
	_b[R:occint4fem] *     occint4fem       +   ///
	_b[R:occint5fem] *     occint5fem       +   ///
	_b[R:occint6fem] *     occint6fem       +	///
	_b[R:occint7fem] *     occint7fem 
}
else if order_expansion==3	{	
	gen a3_interaction=							///
	_b[R:agesq] *          agesq			+	///
	_b[R:expsq] *          expsq		    +	///
	_b[R:agefem] *     agefem		+        	///
	_b[R:agesqfem] *     agesqfem		+    	///
	_b[R:experiencefem] *     experiencefem	+	///
	_b[R:expsqfem] *     expsqfem		+       ///
	_b[R:collegefem] *     collegefem		+   ///
	_b[R:racex2fem] *      racex2fem		+   ///
	_b[R:marriedfem] *     marriedfem		+   ///
	_b[R:widowedfem] *     widowedfem		+   ///
	_b[R:walkra_R_1fem] *  walkra_R_1fem    +	///
	_b[R:dressa_R_1fem] *  dressa_R_1fem    +	///
	_b[R:stoopa_R_1fem] *  stoopa_R_1fem    +	///
	_b[R:beda_R_1fem] *    beda_R_1fem      +	///
	_b[R:bmi_Rfem] *       bmi_Rfem         +   ///
	_b[R:hosp_Rfem] *      hosp_Rfem        +   ///
	_b[R:bsx1fem] *        bsx1fem          +   ///
	_b[R:bsx2fem] *        bsx2fem          +   ///
	_b[R:bsx3fem] *        bsx3fem          +   ///
	_b[R:bsx4fem] *        bsx4fem          +   ///
	_b[R:bsx5fem] *        bsx5fem          +   ///
	_b[R:bsx6fem] *        bsx6fem          +   ///
	_b[R:bsx7fem] *        bsx7fem          +   ///
	_b[R:bsx8fem] *        bsx8fem          +   ///
	_b[R:bsx9fem] *        bsx9fem          +   ///
	_b[R:occint1fem] *     occint1fem       +   ///
	_b[R:occint2fem] *     occint2fem       +   ///
	_b[R:occint3fem] *     occint3fem       +   ///
	_b[R:occint4fem] *     occint4fem       +   ///
	_b[R:occint5fem] *     occint5fem       +   ///
	_b[R:occint6fem] *     occint6fem       +   ///
	_b[R:occint7fem] *     occint7fem +      	///
	_b[R:agecb] *          agecb			+	///
	_b[R:expcb] *          expcb		    +	///
	_b[R:agecbfem] *     agecbfem		+		///
	_b[R:expcbfem] *     expcbfem		
}

gen a3=a3_demo+a3_health+a3_labor+a3_interaction

foreach x in a3 a3_institutional_year a3_demo a3_health a3_labor a3_age	a3_interaction {
	replace `x'=. if R==.
}

nlcom (r21: /c21) (r31: /c31) (r32: /c21*/c31 + sqrt( 1 - /c21^2 )*/c32)
estadd scalar corr_A_L r(b)[1,1]
estadd scalar corr_A_R r(b)[1,2]
estadd scalar corr_L_R r(b)[1,3]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\r(b)[1,1]\r(b)[1,2]\r(b)[1,3]
matrix semoments=semoments\sqrt(r(V)[1,1])\sqrt(r(V)[2,2])\sqrt(r(V)[3,3])
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar rho_aw_ap=-r(b)[1,2]	

gen rho_A_L=r(b)[1,1]
gen rho_A_R=r(b)[1,2]
gen rho_L_R=r(b)[1,3]

preserve
	/*The vignette identifies -X*a(Lbar)=a1_fromV	*/
	/*The work limitations identifies X*a(L*)-X*aL(bar)=a1 (from which need to subtract the female effect) */
	/*To retrieve X*a(L*)=a1+X*aL(bar)=a1-a1_fromV. That's what we have below*/
	gen xbL=(a1-betaFL*female)-(a1_fromV)			 
	gen xbL_only_inst=(a1-betaFL*female)-(a1_inst_fromV)			 
	keep hhidpn a1 a2 xbL* female A R college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R year occ_interm bsx1-bsx9 agesq-occint8fem
	save $data\data_to_est_alpha_nonlin,replace			
	/* Do the saving here, before dropping the first 3 waves for estimating log spending */
restore

qui tab wave,gen(waved) 
qui tab occ_interm,gen(occx)
qui tab year,gen(yeard)


drop rho* bF*

*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global spec_C 	female college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave ///
				i.occ_interm logy lowcash govt_hins priv_hins 		
global spec_U 	$spec_C spin somecov* sp_medexp*		
************************************************************************************************************************************************************************************

gen Aw=1-R

gen A_Aw_3ml=-a3
gen A_App_3ml=a2


*********************************************************
drop if wave<3 			/*no spending data in waves 1,2*/
*********************************************************

duplicates tag hhidpn year logx,gen(check_type12_only)			/*Need 1 obs per year/per HH*/
tab check_type12_only,miss

gen nonzero=x>0
gen uncensored_C=logx!=.

***TABLE 10 + SIGMA_{C}  
eststo bigml3: heckman logx $spec_C, select(nonzero= $spec_U ) nolog
gen logxsample=e(sample)
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[logx:female]\r(table)[1,colsof(r(table))-1]
matrix semoments=semoments\_se[logx:female]\r(table)[2,colsof(r(table))-1]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
scalar sigmac	=r(table)[1,colsof(r(table))-1]
scalar se_sigmac=r(table)[2,colsof(r(table))-1]

predict xb,xb
gen u=logx-xb

drop waved* occx* yeard*

												
biprobit (uncensored_C =$spec_U) (L =	$c_L),nolog	
predict a_unc1,xb1
predict a_L,xb2
scalar rho_uC_uL=e(rho)

gen corr_C_0=normalden(a_unc1)*(normal((1-rho_uC_uL^2)^(-0.5)*(a_L-rho_uC_uL*a_unc1)))*binormal(a_unc1,a_L,rho_uC_uL)^(-1)		/*These expressions are the ones verified via Monte Carlo simulations*/
gen corr_C_L=normalden(a_L)   *(normal((1-rho_uC_uL^2)^(-0.5)*(a_unc1-rho_uC_uL*a_L)))*binormal(a_unc1,a_L,rho_uC_uL)^(-1)

***TABLE 10: RHO_{C,L}
reg u corr_C_0 corr_C_L if L==1 & uncensored_C==1,nocons
gen rcl = _b[corr_C_L]
 
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_L]						/*This is corr(u_C,u_L)*/
matrix semoments=semoments\_se[corr_C_L]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
											
												
biprobit (uncensored_C =$spec_U) (A =	$c_A),nolog	
predict a_unc2,xb1
predict a_App,xb2
scalar rho_uC_uA=e(rho)
scalar rho_unc_ap=e(rho)
	
gen corr_C_0_2=normalden(a_unc2)*(normal((1-rho_uC_uA^2)^(-0.5)*(a_App-rho_uC_uA*a_unc2)))*binormal(a_unc2,a_App,rho_uC_uA)^(-1)	/*These expressions are the ones verified via MC simulations*/
gen corr_C_App=normalden(a_App) *(normal((1-rho_uC_uA^2)^(-0.5)*(a_unc2-rho_uC_uA*a_App)))*binormal(a_unc2,a_App,rho_uC_uA)^(-1)

***TABLE 10: RHO_{C,A}
reg u corr_C_0_2 corr_C_App if A==1 & uncensored_C==1,nocons
gen rca = _b[corr_C_App]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_App]				/*This is corr(u_C,u_A)*/
matrix semoments=semoments\_se[corr_C_App]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar rho_C_0_TBUB		=_b[corr_C_0_2]
scalar rho_C_App_TBUB	=_b[corr_C_App]

*gen Aw=1-R

biprobit (uncensored_C =$spec_U) (Aw =	$c_R),nolog	
predict a_unc3,xb1
predict a_Aw,xb2
scalar rho_unc_aw=e(rho)

egen c2=rmean(a_unc2 a_unc3)
gen  c3=A_App_3ml
gen  c4=A_Aw_3ml
scalar r23=rho_unc_ap
scalar r24=rho_unc_aw 
scalar r34=rho_aw_ap
scalar r34_2=(r34-r23*r24)/(sqrt(1-r23^2)*sqrt(1-r24^2))
scalar r24_3=(r24-r23*r34)/(sqrt(1-r23^2)*sqrt(1-r34^2))
scalar r23_4=(r23-r24*r34)/(sqrt(1-r24^2)*sqrt(1-r34^2))
matrix Sigma = (1, r23, r24 \ r23, 1, r34 \ r24, r34, 1)

mdraws, neq(3) dr(20) prefix(q) random
matrix C=cholesky(Sigma)
egen F3 = mvnp(c2 c3 c4), chol(C) dr(20) prefix(q) adoonly

gen corr_C_00=normalden(c2)*binormal((1-r23^2)^(-0.5)*(c3-r23*c2),(1-r24^2)^(-0.5)*(c4-r24*c2),r34_2)/F3
gen corr_C_ap=normalden(c3)*binormal((1-r23^2)^(-0.5)*(c2-r23*c3),(1-r34^2)^(-0.5)*(c4-r34*c3),r24_3)/F3
gen corr_C_aw=normalden(c4)*binormal((1-r24^2)^(-0.5)*(c2-r24*c4),(1-r34^2)^(-0.5)*(c3-r34*c4),r23_4)/F3

gen v=u-(rho_C_0_TBUB)*corr_C_00-(rho_C_App_TBUB)*corr_C_ap			/*use consistent estimates of corr.coeff. from other parts of the code*/

***TABLE 10: RHO_{C,R}
reg v corr_C_aw if Aw==1 & uncensored_C==1 & A==1,nocons
gen rcr = -_b[corr_C_aw]
scalar b_CA=_b[corr_C_aw]
scalar b_CAse=_se[corr_C_aw]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_aw]
matrix semoments=semoments\_se[corr_C_aw]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

matrix M=moments,semoments
matrix rownames M = 1_bF_V 2_bFv_V 3_corr_L_V 4_b_F_L 5_bF_A 6_corr_L_A 7_corr_R_A 8_corr_L_R 9_bF_C 10_sdu_C 11_EuC_L 12_EuC_A 13_EuC_Aw
matrix M=M[4..6,1..2]\M[1..3,1..2]\M[8..8,1..2]\M[9..11,1..2]\M[13..13,1..2]\M[12..12,1..2]\M[7..7,1..2]

preserve
	matrix list M
	svmat M
	keep M1 M2
	keep if M1!=.
	gen str20 namevar=""
	replace namevar="bF_L"            in 1
	replace namevar="bF_A"            in 2
	replace namevar="corr_L_A"        in 3
	replace namevar="bF_V" 			  in 4
	replace namevar="bFv_V"           in 5
	replace namevar="corr_L_V"        in 6
	replace namevar="corr_L_R"        in 7
	replace namevar="bF_C"            in 8
	replace namevar="sdu_C"           in 9
	replace namevar="EuC_L"           in 10
	replace namevar="EuC_Aw"          in 11
	replace namevar="EuC_A"   		  in 12
	replace namevar="corr_R_A"        in 13
	if order_expansion==2	{	
		save $data\moments_for_R_nonlin_expansion2,replace
	}
	else if order_expansion==3	{	
		save $data\moments_for_R_nonlin_expansion3,replace
	}
	
restore


***To get alpha
cap program drop get_alpha
program define get_alpha
	scalar sigma	=sqrt(1+s_csi^2)			/*This is [1+var(noise_SSA)] */
	gen xD=xbL_only*sqrt(1+s_omega^2+s_psi^2)+mu*female			
	ren a2 xA
	scalar rho_euA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)								/*We want the corr b/w epsilon and uA, not between uL and uA*/
	gen pDi=binormal(xA,xD,rho_euA)/normal(xA)
	gen m1i=normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,xD,rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,xD,rho_euA)^(-1)	
	gen m0i=-normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,-xD,-rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(1-normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,-xD,-rho_euA)^(-1)	
	scalar sigma=sqrt(1+s_csi^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)
	su part1 if female==1 & A==1
	scalar p1_1=r(mean)
	su part1 if female==0 & A==1
	scalar p1_0=r(mean)
	scalar part1m=p1_1-p1_0
	su part2 if female==1 & A==1
	scalar p2_1=r(mean)
	su part2 if female==0 & A==1
	scalar p2_0=r(mean)
	scalar part2m=p2_1-p2_0
	su part3 if female==1 & A==1
	scalar part3m=r(mean)
	gen rho_FX		=part1+part2
	gen lambda_FX	=part3*female
	gen F=female
	probit R female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx1-bsx9 agesq-occint8fem
	predict PI,pr
	gen psi=invnorm(PI)
	su psi 		if F==1 & R!=.
	scalar psi1=r(mean)
	su psi 		if F==0 & R!=.
	scalar psi0=r(mean)
	su rho 		if F==1 & R!=.
	scalar rho1=r(mean)
	su rho 		if F==0 & R!=.
	scalar rho0=r(mean)
	su lambda 	if F==1 & R!=.
	scalar lambda1=r(mean)
	scalar alpha1=((psi1-psi0)-(rho1-rho0))/lambda1		/*New alpha D-in-D: Feb 2024*/
	scalar list alpha1
end

cd $data

u V_pars,clear
sort draw
merge draw using R_pars
drop if _merge!=3
drop _merge
sort draw
merge draw using A_pars
drop if _merge!=3
drop _merge

egen nonconv=rsum(cvg_VL cvg_ALR cvg_C cvg_0L cvg_0A cvg_0Aw_R cvg_partcorr cvg_chol)
keep if nonconv==8

set seed 110968
gen uu=uniform()
sort uu
keep in 1/200

drop uu draw

ren alpha1 alpha_old
ren alpha1_alt alpha1		
su alpha1
scalar se_alpha1=r(sd)

corr bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L E_uC_R E_uC_A rho_A_R alpha1,cov 
matrix V_alpha=r(C)

corr bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L E_uC_R E_uC_A rho_A_R,cov 
matrix V=r(C)
mat M3=vecdiag(V)'


***THIS IS EWMD (with the I matrix as weighting matrix)
if order_expansion==2	{	
	u $data\moments_for_R_nonlin_expansion2,clear
}
else if order_expansion==3	{	
	u $data\moments_for_R_nonlin_expansion3,clear
}
mkmat M1
mata: mata clear
mata:
void eval0(todo, b, v, g, H)
 {
		  x1 =(b[4]-b[5])/sqrt(1+b[1]^2+b[9]^2)										
		  x2 =(b[4]-b[5]-b[6])/sqrt(1+b[1]^2+b[9]^2+b[2]^2)							
		  x3 =(1+b[1]^2+b[9]^2)/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2+b[2]^2))	
		  x4 =-(b[5]-b[7])/sqrt(1+b[1]^2+b[9]^2)											
		  x5 =(b[4])/sqrt(1+b[1]^2+b[9]^2)											
		  x6 =b[9]^2/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2))					
		  x7 =-1/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[11]^2))								
		  x8 =b[8]*b[4]+b[10]															
		  x9 =sqrt(b[8]^2+b[3]^2)														
		  x10=b[8]/sqrt(1+b[1]^2+b[9]^2)												
		  x11=b[8]/sqrt(1+b[11]^2)
		  x12=b[8]/sqrt(1+b[1]^2+b[9]^2+b[2]^2)
		  x13=-1/(sqrt(1+b[1]^2+b[9]^2+b[2]^2)*sqrt(1+b[11]^2))
		  M1=st_matrix("M1")
		  M2=st_matrix("M3")
		  m1 =M1[1,1]
		  m2 =M1[2,1] 
		  m3 =M1[3,1]
		  m4 =M1[4,1]
		  m5 =M1[5,1] 
		  m6 =M1[6,1]
		  m7 =M1[7,1] 
		  m8 =M1[8,1]
		  m9 =M1[9,1]
		  m10=M1[10,1]
		  m11=M1[11,1]
		  m12=M1[12,1]
		  m13=M1[13,1]
		  v1 =M2[1,1]
		  v2 =M2[2,1] 
		  v3 =M2[3,1]
		  v4 =M2[4,1]
		  v5 =M2[5,1]
		  v6 =M2[6,1]
		  v7 =M2[7,1] 
		  v8 =M2[8,1]
		  v9 =M2[9,1]
		  v10=M2[10,1]
		  v11=M2[11,1]
		  v12=M2[12,1]
		  v13=M2[13,1]
	v=((x1-m1)^2/1)+((x2-m2)^2/1)+((x3-m3)^2/1)+((x4-m4)^2/1)+((x5-m5)^2/1)+((x6-m6)^2/1)+((x7-m7)^2/1)+((x8-m8)^2/1)+((x9-m9)^2/1)+((x10-m10)^2/1)+((x11-m11)^2/1)+((x12-m12)^2/1)+((x13-m13)^2/1)
 }
S = optimize_init()
optimize_init_evaluator(S, &eval0())
optimize_init_evaluatortype(S, "v0")
optimize_init_params(S, (1,1,1,1,1,1,1,1,1,1,1))
optimize_init_which(S, "min")
optimize_init_which(S, "min")
optimize_init_conv_ptol(S, 1e-16)
optimize_init_conv_vtol(S, 1e-16)
b = optimize(S)
b'
st_matrix("beta_ewmd",b')
end


scalar mu		=beta_ewmd[4,1]
scalar s_omega	=beta_ewmd[1,1]
scalar s_psi	=beta_ewmd[9,1]	
scalar s_v		=beta_ewmd[2,1]
scalar s_csi	=beta_ewmd[11,1]	
u data_to_est_alpha_nonlin,clear
get_alpha
matrix beta=beta_ewmd\alpha1
svmat beta
keep beta
drop if beta==.
ren beta1 beta
gen str20 namepar=""
replace namepar="sigma_omega"      in 1
replace namepar="sigma_v"   	   in 2
replace namepar="sigma_phi"        in 3
replace namepar="pi"  		       in 4
replace namepar="gamma"  	       in 5
replace namepar="tau"   	       in 6
replace namepar="theta"  	       in 7
replace namepar="lamda"  	       in 8
replace namepar="sigma_psi"        in 9
replace namepar="psi"   	       in 10
replace namepar="sigma_csi"        in 11
replace namepar="alpha"        	   in 12
order namepar beta

*****************************
if order_expansion==2	{	
	save results_nonlin_expansion2,replace
}
else if order_expansion==3	{	
	save results_nonlin_expansion3,replace
}


***Now Run this for order_expansion=2*/
********************************************************************************
scalar define order_expansion=2		/*2=quadratic expansion, 3=cubic expansion*/
********************************************************************************

cd $data

************************************************************************************************************************************************************************************
******VIGNETTE DATA
************************************************************************************************************************************************************************************
************************************************************************************************************************************************************************************
u $data\vignette_data, clear
******************************************************************************************************************************************************************

ren walkra walkra_R_1 
ren dressa dressa_R_1 
ren stoopa stoopa_R_1 
ren beda   beda_R_1   
ren bmi    bmi_R      
ren hosp   hosp_R      

qui tab occ_interm,gen(occx)
qui tab vign_id,gen(vigndomx)

************************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global spec_L 		female college racex2 doctor* married widowed age  walkra dressa stoopa beda hosp bmi i.occ_interm experience
global spec_V 		$spec_L vign_f i.vign_id 
************************************************************************************************************************************************************************************

******************************************************************************************************************************************************************
**Sample selection (drop those w/ missing data in the variables used for the 2 regressions)
qui biprobit (vign_severe = $spec_V) (cantwork = $spec_L) if year==2007, vce(cluster hhidpn)	
keep if e(sample)
******************************************************************************************************************************************************************

*This is the bivariate probit L(respondent), L(vignette)
****TABLE 9 + \RHO_{L,V}
eststo bigml2: biprobit (vign_severe = $spec_V) (cantwork = $spec_L) if year==2007, vce(cluster hhidpn)	

matrix bV=e(b)
matrix bV=bV'

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments	=_b[vign_severe:female]\_b[vign_severe:vign_f]\r(table)[1,colsof(r(table))]
matrix semoments=_se[vign_severe:female]\_se[vign_severe:vign_f]\r(table)[2,colsof(r(table))]	
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
************************************************************************************************************************************************

************************************************************************************************************************************************
estadd scalar corr_L_V 	r(table)[1,colsof(r(table))]


************************************************************************************************************************************************************************************
********REJECTION DATA *************************************************************************************************************************************************************************************
*************************************************************************************************************************************************************************************

u $data\rejection_data, clear
***QUADRATIC EXPANSION
g agesq=age^2
g expsq=experience^2

foreach x of varlist age experience agesq expsq college racex2 married widowed walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R {
	gen `x'fem=`x'*female
}

*****************************
if order_expansion==3	{	
	gen agecb=age^3
	gen expcb=experience^3
	gen agecbfem=agecb*female
	gen expcbfem=expcb*female
}
*****************************

forvalues j=1(1)9 {
	gen bsx`j'fem=bsx`j'*female
}
forvalues j=1(1)8 {
	gen occint`j'fem=(occ_interm==`j')*female
}
save $data\rejection_data_extended,replace

u $data\rejection_data_extended, clear

/*This is -X*a(Lbar)*/
#delimit ;
gen a1_fromV=
bV[2 ,1]		*college+  
bV[3 ,1]         *racex2  +
bV[4 ,1]   *doctor_told_has_hbp +
bV[5 ,1]   *doctor_told_has_psy +
bV[6 ,1]   *doctor_told_has_hea +
bV[7 ,1]   *doctor_told_has_art +
bV[8 ,1]   *doctor_told_has_dia +
bV[9 ,1]   *doctor_told_has_lun +
bV[10,1]   *doctor_told_has_str +
bV[11,1]   *doctor_told_has_can +
bV[12,1]         *married +
bV[13,1]         *widowed +
bV[14,1]             *age +
bV[15,1]      *walkra_R_1 +
bV[16,1]      *dressa_R_1 +
bV[17,1]      *stoopa_R_1 +
bV[18,1]        *beda_R_1 +
bV[19,1]          *hosp_R +
bV[20,1]          *bmi_R  +
bV[21,1]     *(occ_interm==1)+
bV[22,1]     *(occ_interm==2)+
bV[23,1]     *(occ_interm==3)+
bV[24,1]     *(occ_interm==4)+
bV[25,1]     *(occ_interm==5)+
bV[26,1]     *(occ_interm==6)+
bV[27,1]     *(occ_interm==7)+
bV[28,1]     *(occ_interm==8)+
bV[29,1]      *experience;
#delimit cr

/*This is -X*a(Lbar)(institutional variables only)*/
#delimit ;
gen a1_inst_fromV=
bV[2 ,1]		*college+  
bV[4 ,1]   *doctor_told_has_hbp +
bV[5 ,1]   *doctor_told_has_psy +
bV[6 ,1]   *doctor_told_has_hea +
bV[7 ,1]   *doctor_told_has_art +
bV[8 ,1]   *doctor_told_has_dia +
bV[9 ,1]   *doctor_told_has_lun +
bV[10,1]   *doctor_told_has_str +
bV[11,1]   *doctor_told_has_can +
bV[14,1]             *age +
bV[15,1]      *walkra_R_1 +
bV[16,1]      *dressa_R_1 +
bV[17,1]      *stoopa_R_1 +
bV[18,1]        *beda_R_1 +
bV[19,1]          *hosp_R +
bV[20,1]          *bmi_R  +bV[21,1]     *(occ_interm==1)+
bV[22,1]     *(occ_interm==2)+
bV[23,1]     *(occ_interm==3)+
bV[24,1]     *(occ_interm==4)+
bV[25,1]     *(occ_interm==5)+
bV[26,1]     *(occ_interm==6)+
bV[27,1]     *(occ_interm==7)+
bV[28,1]     *(occ_interm==8)+
bV[29,1]      *experience;
#delimit cr


************************************************************************************************************************************************************************************
*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global c_L 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm 
global c_A 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave i.occ_interm logbs 
global c_R 	female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx1-bsx9 agesq-occint8fem 
************************************************************************************************************************************************************************************

gen A=applied
gen R=nosuccess
gen L=cantwork_notemp

gen 	index2=1 if R==1 & A==1
replace index2=2 if R==0 & A==1
replace index2=3 if A==0

set seed 2134569

// Create Halton draw variables
mdraws, dr(10) neq(3) prefix(z)			/*If dr(N), then in egen `sp3' and egen `sp2' adjust dr(N) accordingly*/

// Get initial values
quietly {
        probit A $c_A
        mat b1 = e(b)
        mat coleq b1 = A
        probit L $c_L
        mat b2 = e(b)
        mat coleq b2 = L
        probit R $c_R
        mat b3 = e(b)
        mat coleq b3 = R
        mat b0 = b1, b2, b3
}


cap program drop myll2
program define myll2
args lnf xb1 xb2 xb3 c21 c31 c32
        tempvar sp2 sp3 k1 k2 k3
        quietly {
            gen double `k1' = 2*$ML_y1 - 1
            gen double `k2' = 2*$ML_y2 - 1
            gen double `k3' = 2*$ML_y3 - 1
            tempname cf21 cf31 cf32 cf22 cf33 C1 C2
            su `c21', meanonly
            scalar `cf21' = r(mean)
            su `c31', meanonly
            scalar `cf31' = r(mean)
            su `c32', meanonly
            scalar `cf32' = r(mean)
            scalar `cf22' = sqrt( 1 - `cf21'^2 )
            scalar `cf33' = sqrt( 1 - `cf31'^2 - `cf32'^2 )
            tempname C1 C2
            mat `C1' = (1, 0, 0 \ `cf21', `cf22', 0 \ `cf31', `cf32', `cf33')
            mat `C2' = (1, 0 \ `cf21', `cf22' )
            egen `sp3' = mvnp(`xb1' `xb2' `xb3') if $ML_y1==1, ///
                    chol(`C1') dr(10) prefix(z) signs(`k1' `k2' `k3') adoonly
            egen `sp2' = mvnp(`xb1' `xb2' ) if $ML_y1==0, ///
                    chol(`C2') dr(10) prefix(z) signs(`k1' `k2') adoonly
            replace `lnf'= ln(`sp3') if $ML_y1==1
            replace `lnf'= ln(`sp2') if $ML_y1==0
        }
end

ml model lf myll2 	(A: A=$c_A) ///
					(L: L=$c_L) ///
					(R: R=$c_R) ///
					/c21 /c31 /c32 ///
					, missing

ml init b0
eststo bigml1:ml maximize

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[L:female]\_b[A:female]
matrix semoments=semoments\_se[L:female]\_se[A:female]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar betaFL=_b[L:female]	
scalar delta =_b[R:female]	/*to be used for pinning down \alpha*/

gen bF_L=_b[L:female]
gen bF_A=_b[A:female] 
gen bF_R=_b[R:female]
	
#delimit;
gen a1_age=_b[L:age]*age;
gen a1_institutional_year=
_b[L:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[L:doctor_told_has_psy]* doctor_told_has_psy+
_b[L:doctor_told_has_hea]* doctor_told_has_hea+
_b[L:doctor_told_has_art]* doctor_told_has_art+
_b[L:doctor_told_has_dia]* doctor_told_has_dia+
_b[L:doctor_told_has_lun]* doctor_told_has_lun+
_b[L:doctor_told_has_str]* doctor_told_has_str+
_b[L:doctor_told_has_can]* doctor_told_has_can+
_b[L:walkra_R_1]*             walkra_R_1+
_b[L:dressa_R_1]*             dressa_R_1+
_b[L:stoopa_R_1]*             stoopa_R_1+
_b[L:beda_R_1]*                 beda_R_1+
_b[L:bmi_R]*                       bmi_R+
_b[L:hosp_R]*                     hosp_R+
_b[L:age]*                           age+
_b[L:college]*                   college+
_b[L:experience]*             experience+
_b[L:1b.occ_interm]*        (occ_interm==1)+
_b[L:2.occ_interm]*         (occ_interm==2)+
_b[L:3.occ_interm]*         (occ_interm==3)+
_b[L:4.occ_interm]*         (occ_interm==4)+
_b[L:5.occ_interm]*         (occ_interm==5)+
_b[L:6.occ_interm]*         (occ_interm==6)+
_b[L:7.occ_interm]*         (occ_interm==7)+
_b[L:8.occ_interm]*         (occ_interm==8)+
_b[L:1b.wave]*                  (wave==1)+
_b[L:2.wave]*                   (wave==2)+
_b[L:3.wave]*                   (wave==3)+
_b[L:4.wave]*                   (wave==4)+
_b[L:5.wave]*                   (wave==5)+
_b[L:6.wave]*                   (wave==6)+
_b[L:7.wave]*                   (wave==7)+
_b[L:8.wave]*                   (wave==8)+
_b[L:9.wave]*                   (wave==9)+
_b[L:10.wave]*                  (wave==10)+
_b[L:11.wave]*                  (wave==11)+
_b[L:12.wave]*                  (wave==12)+
_b[L:13.wave]*                  (wave==13)+
_b[L:14.wave]*                  (wave==14)+
_b[L:15.wave]*                  (wave==15);
gen a1_health=
_b[L:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[L:doctor_told_has_psy]* doctor_told_has_psy+
_b[L:doctor_told_has_hea]* doctor_told_has_hea+
_b[L:doctor_told_has_art]* doctor_told_has_art+
_b[L:doctor_told_has_dia]* doctor_told_has_dia+
_b[L:doctor_told_has_lun]* doctor_told_has_lun+
_b[L:doctor_told_has_str]* doctor_told_has_str+
_b[L:doctor_told_has_can]* doctor_told_has_can+
_b[L:walkra_R_1]*             walkra_R_1+
_b[L:dressa_R_1]*             dressa_R_1+
_b[L:stoopa_R_1]*             stoopa_R_1+
_b[L:beda_R_1]*                 beda_R_1+
_b[L:bmi_R]*                       bmi_R+
_b[L:hosp_R]*                     hosp_R;
gen a1_demo=
_b[L:female]*				      female+
_b[L:age]*                           age+
_b[L:racex2]*                     racex2+
_b[L:married]*                   married+
_b[L:widowed]*                   widowed+
_b[L:_cons];
gen a1_labor=
_b[L:college]*                   college+
_b[L:experience]*             experience+
_b[L:1b.wave]*                  (wave==1)+
_b[L:2.wave]*                   (wave==2)+
_b[L:3.wave]*                   (wave==3)+
_b[L:4.wave]*                   (wave==4)+
_b[L:5.wave]*                   (wave==5)+
_b[L:6.wave]*                   (wave==6)+
_b[L:7.wave]*                   (wave==7)+
_b[L:8.wave]*                   (wave==8)+
_b[L:9.wave]*                   (wave==9)+
_b[L:10.wave]*                  (wave==10)+
_b[L:11.wave]*                  (wave==11)+
_b[L:12.wave]*                  (wave==12)+
_b[L:13.wave]*                  (wave==13)+
_b[L:14.wave]*                  (wave==14)+
_b[L:15.wave]*                  (wave==15)+
_b[L:1b.occ_interm]*        (occ_interm==1)+
_b[L:2.occ_interm]*         (occ_interm==2)+
_b[L:3.occ_interm]*         (occ_interm==3)+
_b[L:4.occ_interm]*         (occ_interm==4)+
_b[L:5.occ_interm]*         (occ_interm==5)+
_b[L:6.occ_interm]*         (occ_interm==6)+
_b[L:7.occ_interm]*         (occ_interm==7)+
_b[L:8.occ_interm]*         (occ_interm==8);
#delimit cr
gen 	a1=a1_demo+a1_health+a1_labor
foreach x in a1 a1_institutional_year a1_demo a1_health a1_labor a1_age	{
	replace `x'=. if L==.
}

#delimit;
gen a2_age=_b[A:age]*age;
gen a2_institutional_year=
_b[A:age]*                           age+
_b[A:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[A:doctor_told_has_psy]* doctor_told_has_psy+
_b[A:doctor_told_has_hea]* doctor_told_has_hea+
_b[A:doctor_told_has_art]* doctor_told_has_art+
_b[A:doctor_told_has_dia]* doctor_told_has_dia+
_b[A:doctor_told_has_lun]* doctor_told_has_lun+
_b[A:doctor_told_has_str]* doctor_told_has_str+
_b[A:doctor_told_has_can]* doctor_told_has_can+
_b[A:walkra_R_1]*             walkra_R_1+
_b[A:dressa_R_1]*             dressa_R_1+
_b[A:stoopa_R_1]*             stoopa_R_1+
_b[A:beda_R_1]*                 beda_R_1+
_b[A:bmi_R]*                       bmi_R+
_b[A:hosp_R]*                     hosp_R+
_b[A:college]*                   college+
_b[A:experience]*             experience+
_b[A:1b.occ_interm]*        (occ_interm==1)+
_b[A:2.occ_interm]*         (occ_interm==2)+
_b[A:3.occ_interm]*         (occ_interm==3)+
_b[A:4.occ_interm]*         (occ_interm==4)+
_b[A:5.occ_interm]*         (occ_interm==5)+
_b[A:6.occ_interm]*         (occ_interm==6)+
_b[A:7.occ_interm]*         (occ_interm==7)+
_b[A:8.occ_interm]*         (occ_interm==8)+
_b[A:1b.wave]*                  (wave==1)+
_b[A:2.wave]*                   (wave==2)+
_b[A:3.wave]*                   (wave==3)+
_b[A:4.wave]*                   (wave==4)+
_b[A:5.wave]*                   (wave==5)+
_b[A:6.wave]*                   (wave==6)+
_b[A:7.wave]*                   (wave==7)+
_b[A:8.wave]*                   (wave==8)+
_b[A:9.wave]*                   (wave==9)+
_b[A:10.wave]*                  (wave==10)+
_b[A:11.wave]*                  (wave==11)+
_b[A:12.wave]*                  (wave==12)+
_b[A:13.wave]*                  (wave==13)+
_b[A:14.wave]*                  (wave==14)+
_b[A:15.wave]*                  (wave==15);
gen a2_demo=
_b[A:female]*				      female+
_b[A:age]*                           age+
_b[A:racex2]*                     racex2+
_b[A:married]*                   married+
_b[A:widowed]*                   widowed+
_b[A:logbs]*logbs						+
_b[A:_cons];
gen a2_health=
_b[A:doctor_told_has_hbp]* doctor_told_has_hbp+
_b[A:doctor_told_has_psy]* doctor_told_has_psy+
_b[A:doctor_told_has_hea]* doctor_told_has_hea+
_b[A:doctor_told_has_art]* doctor_told_has_art+
_b[A:doctor_told_has_dia]* doctor_told_has_dia+
_b[A:doctor_told_has_lun]* doctor_told_has_lun+
_b[A:doctor_told_has_str]* doctor_told_has_str+
_b[A:doctor_told_has_can]* doctor_told_has_can+
_b[A:walkra_R_1]*             walkra_R_1+
_b[A:dressa_R_1]*             dressa_R_1+
_b[A:stoopa_R_1]*             stoopa_R_1+
_b[A:beda_R_1]*                 beda_R_1+
_b[A:bmi_R]*                       bmi_R+
_b[A:hosp_R]*                     hosp_R;
gen a2_labor=
_b[A:college]*                   college+
_b[A:experience]*             experience+
_b[A:1b.wave]*                  (wave==1)+
_b[A:2.wave]*                   (wave==2)+
_b[A:3.wave]*                   (wave==3)+
_b[A:4.wave]*                   (wave==4)+
_b[A:5.wave]*                   (wave==5)+
_b[A:6.wave]*                   (wave==6)+
_b[A:7.wave]*                   (wave==7)+
_b[A:8.wave]*                   (wave==8)+
_b[A:9.wave]*                   (wave==9)+
_b[A:10.wave]*                  (wave==10)+
_b[A:11.wave]*                  (wave==11)+
_b[A:12.wave]*                  (wave==12)+
_b[A:13.wave]*                  (wave==13)+
_b[A:14.wave]*                  (wave==14)+
_b[A:15.wave]*                  (wave==15)+
_b[A:1b.occ_interm]*        (occ_interm==1)+
_b[A:2.occ_interm]*         (occ_interm==2)+
_b[A:3.occ_interm]*         (occ_interm==3)+
_b[A:4.occ_interm]*         (occ_interm==4)+
_b[A:5.occ_interm]*         (occ_interm==5)+
_b[A:6.occ_interm]*         (occ_interm==6)+
_b[A:7.occ_interm]*         (occ_interm==7)+
_b[A:8.occ_interm]*         (occ_interm==8);
#delimit cr
gen 	a2=a2_demo+a2_health+a2_labor
foreach x in a2 a2_institutional_year a2_demo a2_health a2_labor a2_age	{
	replace `x'=. if A==.
}

#delimit;
gen a3_age=_b[R:age]*age;
gen a3_institutional_year=
	_b[R:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[R:doctor_told_has_psy]* doctor_told_has_psy+
	_b[R:doctor_told_has_hea]* doctor_told_has_hea+
	_b[R:doctor_told_has_art]* doctor_told_has_art+
	_b[R:doctor_told_has_dia]* doctor_told_has_dia+
	_b[R:doctor_told_has_lun]* doctor_told_has_lun+
	_b[R:doctor_told_has_str]* doctor_told_has_str+
	_b[R:doctor_told_has_can]* doctor_told_has_can+
	_b[R:walkra_R_1]*             walkra_R_1+
	_b[R:dressa_R_1]*             dressa_R_1+
	_b[R:stoopa_R_1]*             stoopa_R_1+
	_b[R:beda_R_1]*                 beda_R_1+
	_b[R:bmi_R]*                       bmi_R+
	_b[R:hosp_R]*                     hosp_R+
	_b[R:bsx1]*bsx1+
	_b[R:bsx2]*bsx2+
	_b[R:bsx3]*bsx3+
	_b[R:bsx4]*bsx4+
	_b[R:bsx5]*bsx5+
	_b[R:bsx6]*bsx6+
	_b[R:bsx7]*bsx7+
	_b[R:bsx8]*bsx8+
	_b[R:bsx9]*bsx9+
	_b[R:age]*                           age+
	_b[R:college]*                   college+
	_b[R:experience]*             experience+
	_b[R:1b.occ_interm]*        (occ_interm==1)+ 
	_b[R:2.occ_interm]*         (occ_interm==2)+
	_b[R:3.occ_interm]*         (occ_interm==3)+
	_b[R:4.occ_interm]*         (occ_interm==4)+
	_b[R:5.occ_interm]*         (occ_interm==5)+
	_b[R:6.occ_interm]*         (occ_interm==6)+
	_b[R:7.occ_interm]*         (occ_interm==7)+
	_b[R:8.occ_interm]*         (occ_interm==8)+
	_b[R:1991b.year]*(year==1991)+
	_b[R:1992.year] *(year==1992)+
	_b[R:1993.year] *(year==1993)+
	_b[R:1994.year] *(year==1994)+
	_b[R:1995.year] *(year==1995)+
	_b[R:1996.year] *(year==1996)+
	_b[R:1997.year] *(year==1997)+
	_b[R:1998.year] *(year==1998)+
	_b[R:1999.year] *(year==1999)+
	_b[R:2000.year] *(year==2000)+
	_b[R:2001.year] *(year==2001)+
	_b[R:2002.year] *(year==2002)+
	_b[R:2003.year] *(year==2003)+
	_b[R:2004.year] *(year==2004)+
	_b[R:2005.year] *(year==2005)+
	_b[R:2006.year] *(year==2006)+
	_b[R:2007.year] *(year==2007)+
	_b[R:2008.year] *(year==2008)+
	_b[R:2009.year] *(year==2009)+
	_b[R:2010.year] *(year==2010)+
	_b[R:2011.year] *(year==2011)+
	_b[R:2012.year] *(year==2012)+
	_b[R:2013.year] *(year==2013)+
	_b[R:2014.year] *(year==2014)+
	_b[R:2015.year] *(year==2015)+
	_b[R:2016.year] *(year==2016)+
	_b[R:2017.year] *(year==2017)+
	_b[R:2018.year] *(year==2018)+
	_b[R:2019.year] *(year==2019)+
	_b[R:2020.year] *(year==2020);	
gen a3_health=
	_b[R:doctor_told_has_hbp]* doctor_told_has_hbp+
	_b[R:doctor_told_has_psy]* doctor_told_has_psy+
	_b[R:doctor_told_has_hea]* doctor_told_has_hea+
	_b[R:doctor_told_has_art]* doctor_told_has_art+
	_b[R:doctor_told_has_dia]* doctor_told_has_dia+
	_b[R:doctor_told_has_lun]* doctor_told_has_lun+
	_b[R:doctor_told_has_str]* doctor_told_has_str+
	_b[R:doctor_told_has_can]* doctor_told_has_can+
	_b[R:walkra_R_1]*             walkra_R_1+
	_b[R:dressa_R_1]*             dressa_R_1+
	_b[R:stoopa_R_1]*             stoopa_R_1+
	_b[R:beda_R_1]*                 beda_R_1+
	_b[R:bmi_R]*                       bmi_R+
	_b[R:hosp_R]*                     hosp_R+
	_b[R:bsx1]*bsx1+
	_b[R:bsx2]*bsx2+
	_b[R:bsx3]*bsx3+
	_b[R:bsx4]*bsx4+
	_b[R:bsx5]*bsx5+
	_b[R:bsx6]*bsx6+
	_b[R:bsx7]*bsx7+
	_b[R:bsx8]*bsx8+
	_b[R:bsx9]*bsx9;
gen a3_demo=
	_b[R:female]*				      female+
	_b[R:age]*                           age+
	_b[R:racex2]*                     racex2+
	_b[R:married]*                   married+
	_b[R:widowed]*                   widowed+
	_b[R:_cons];
gen a3_labor=
	_b[R:college]*                   college+
	_b[R:experience]*             experience+
	_b[R:1991b.year]*(year==1991)+
	_b[R:1992.year] *(year==1992)+
	_b[R:1993.year] *(year==1993)+
	_b[R:1994.year] *(year==1994)+
	_b[R:1995.year] *(year==1995)+
	_b[R:1996.year] *(year==1996)+
	_b[R:1997.year] *(year==1997)+
	_b[R:1998.year] *(year==1998)+
	_b[R:1999.year] *(year==1999)+
	_b[R:2000.year] *(year==2000)+
	_b[R:2001.year] *(year==2001)+
	_b[R:2002.year] *(year==2002)+
	_b[R:2003.year] *(year==2003)+
	_b[R:2004.year] *(year==2004)+
	_b[R:2005.year] *(year==2005)+
	_b[R:2006.year] *(year==2006)+
	_b[R:2007.year] *(year==2007)+
	_b[R:2008.year] *(year==2008)+
	_b[R:2009.year] *(year==2009)+
	_b[R:2010.year] *(year==2010)+
	_b[R:2011.year] *(year==2011)+
	_b[R:2012.year] *(year==2012)+
	_b[R:2013.year] *(year==2013)+
	_b[R:2014.year] *(year==2014)+
	_b[R:2015.year] *(year==2015)+
	_b[R:2016.year] *(year==2016)+
	_b[R:2017.year] *(year==2017)+
	_b[R:2018.year] *(year==2018)+
	_b[R:2019.year] *(year==2019)+
	_b[R:2020.year] *(year==2020)+
	_b[R:1b.occ_interm]*        (occ_interm==1)+ 
	_b[R:2.occ_interm]*         (occ_interm==2)+
	_b[R:3.occ_interm]*         (occ_interm==3)+
	_b[R:4.occ_interm]*         (occ_interm==4)+
	_b[R:5.occ_interm]*         (occ_interm==5)+
	_b[R:6.occ_interm]*         (occ_interm==6)+
	_b[R:7.occ_interm]*         (occ_interm==7)+
	_b[R:8.occ_interm]*         (occ_interm==8);
#delimit cr

if order_expansion==2	{	
	gen a3_interaction= 					 	///
	_b[R:agesq] *          agesq			+	///	
	_b[R:expsq] *          expsq		    +	///
	_b[R:agefem] *     agefem		+			///
	_b[R:agesqfem] *     agesqfem		+		///
	_b[R:experiencefem] *     experiencefem	+   ///
	_b[R:expsqfem] *     expsqfem		+       ///
	_b[R:collegefem] *     collegefem		+   ///
	_b[R:racex2fem] *      racex2fem		+   ///
	_b[R:marriedfem] *     marriedfem		+   ///
	_b[R:widowedfem] *     widowedfem		+	///
	_b[R:walkra_R_1fem] *  walkra_R_1fem    +	///
	_b[R:dressa_R_1fem] *  dressa_R_1fem    +	///
	_b[R:stoopa_R_1fem] *  stoopa_R_1fem    +	///
	_b[R:beda_R_1fem] *    beda_R_1fem      +   ///
	_b[R:bmi_Rfem] *       bmi_Rfem         +   ///
	_b[R:hosp_Rfem] *      hosp_Rfem        +   ///
	_b[R:bsx1fem] *        bsx1fem          +   ///
	_b[R:bsx2fem] *        bsx2fem          +   ///
	_b[R:bsx3fem] *        bsx3fem          +   ///
	_b[R:bsx4fem] *        bsx4fem          +   ///
	_b[R:bsx5fem] *        bsx5fem          +   ///
	_b[R:bsx6fem] *        bsx6fem          +   ///
	_b[R:bsx7fem] *        bsx7fem          +   ///
	_b[R:bsx8fem] *        bsx8fem          +   ///
	_b[R:bsx9fem] *        bsx9fem          +   ///
	_b[R:occint1fem] *     occint1fem       +   ///
	_b[R:occint2fem] *     occint2fem       +   ///
	_b[R:occint3fem] *     occint3fem       +   ///
	_b[R:occint4fem] *     occint4fem       +   ///
	_b[R:occint5fem] *     occint5fem       +   ///
	_b[R:occint6fem] *     occint6fem       +	///
	_b[R:occint7fem] *     occint7fem 
}
else if order_expansion==3	{	
	gen a3_interaction=							///
	_b[R:agesq] *          agesq			+	///
	_b[R:expsq] *          expsq		    +	///
	_b[R:agefem] *     agefem		+        	///
	_b[R:agesqfem] *     agesqfem		+    	///
	_b[R:experiencefem] *     experiencefem	+	///
	_b[R:expsqfem] *     expsqfem		+       ///
	_b[R:collegefem] *     collegefem		+   ///
	_b[R:racex2fem] *      racex2fem		+   ///
	_b[R:marriedfem] *     marriedfem		+   ///
	_b[R:widowedfem] *     widowedfem		+   ///
	_b[R:walkra_R_1fem] *  walkra_R_1fem    +	///
	_b[R:dressa_R_1fem] *  dressa_R_1fem    +	///
	_b[R:stoopa_R_1fem] *  stoopa_R_1fem    +	///
	_b[R:beda_R_1fem] *    beda_R_1fem      +	///
	_b[R:bmi_Rfem] *       bmi_Rfem         +   ///
	_b[R:hosp_Rfem] *      hosp_Rfem        +   ///
	_b[R:bsx1fem] *        bsx1fem          +   ///
	_b[R:bsx2fem] *        bsx2fem          +   ///
	_b[R:bsx3fem] *        bsx3fem          +   ///
	_b[R:bsx4fem] *        bsx4fem          +   ///
	_b[R:bsx5fem] *        bsx5fem          +   ///
	_b[R:bsx6fem] *        bsx6fem          +   ///
	_b[R:bsx7fem] *        bsx7fem          +   ///
	_b[R:bsx8fem] *        bsx8fem          +   ///
	_b[R:bsx9fem] *        bsx9fem          +   ///
	_b[R:occint1fem] *     occint1fem       +   ///
	_b[R:occint2fem] *     occint2fem       +   ///
	_b[R:occint3fem] *     occint3fem       +   ///
	_b[R:occint4fem] *     occint4fem       +   ///
	_b[R:occint5fem] *     occint5fem       +   ///
	_b[R:occint6fem] *     occint6fem       +   ///
	_b[R:occint7fem] *     occint7fem +      	///
	_b[R:agecb] *          agecb			+	///
	_b[R:expcb] *          expcb		    +	///
	_b[R:agecbfem] *     agecbfem		+		///
	_b[R:expcbfem] *     expcbfem		
}

gen a3=a3_demo+a3_health+a3_labor+a3_interaction

foreach x in a3 a3_institutional_year a3_demo a3_health a3_labor a3_age	a3_interaction {
	replace `x'=. if R==.
}

nlcom (r21: /c21) (r31: /c31) (r32: /c21*/c31 + sqrt( 1 - /c21^2 )*/c32)
estadd scalar corr_A_L r(b)[1,1]
estadd scalar corr_A_R r(b)[1,2]
estadd scalar corr_L_R r(b)[1,3]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\r(b)[1,1]\r(b)[1,2]\r(b)[1,3]
matrix semoments=semoments\sqrt(r(V)[1,1])\sqrt(r(V)[2,2])\sqrt(r(V)[3,3])
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar rho_aw_ap=-r(b)[1,2]	

gen rho_A_L=r(b)[1,1]
gen rho_A_R=r(b)[1,2]
gen rho_L_R=r(b)[1,3]

preserve
	/*The vignette identifies -X*a(Lbar)=a1_fromV	*/
	/*The work limitations identifies X*a(L*)-X*aL(bar)=a1 (from which need to subtract the female effect) */
	/*To retrieve X*a(L*)=a1+X*aL(bar)=a1-a1_fromV. That's what we have below*/
	gen xbL=(a1-betaFL*female)-(a1_fromV)			 
	gen xbL_only_inst=(a1-betaFL*female)-(a1_inst_fromV)			 
	keep hhidpn a1 a2 xbL* female A R college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R year occ_interm bsx1-bsx9 agesq-occint8fem
	save $data\data_to_est_alpha_nonlin,replace			
	/* Do the saving here, before dropping the first 3 waves for estimating log spending */
restore

qui tab wave,gen(waved) 
qui tab occ_interm,gen(occx)
qui tab year,gen(yeard)


drop rho* bF*

*****************************EMPIRICAL SPECIFICATIONS ******************************************************************************************************************************
global spec_C 	female college racex2 age experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.wave ///
				i.occ_interm logy lowcash govt_hins priv_hins 		
global spec_U 	$spec_C spin somecov* sp_medexp*		
************************************************************************************************************************************************************************************

gen Aw=1-R

gen A_Aw_3ml=-a3
gen A_App_3ml=a2


*********************************************************
drop if wave<3 			/*no spending data in waves 1,2*/
*********************************************************

duplicates tag hhidpn year logx,gen(check_type12_only)			/*Need 1 obs per year/per HH*/
tab check_type12_only,miss

gen nonzero=x>0
gen uncensored_C=logx!=.

***TABLE 10 + SIGMA_{C}  
eststo bigml3: heckman logx $spec_C, select(nonzero= $spec_U ) nolog
gen logxsample=e(sample)
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[logx:female]\r(table)[1,colsof(r(table))-1]
matrix semoments=semoments\_se[logx:female]\r(table)[2,colsof(r(table))-1]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
scalar sigmac	=r(table)[1,colsof(r(table))-1]
scalar se_sigmac=r(table)[2,colsof(r(table))-1]

predict xb,xb
gen u=logx-xb

drop waved* occx* yeard*

												
biprobit (uncensored_C =$spec_U) (L =	$c_L),nolog	
predict a_unc1,xb1
predict a_L,xb2
scalar rho_uC_uL=e(rho)

gen corr_C_0=normalden(a_unc1)*(normal((1-rho_uC_uL^2)^(-0.5)*(a_L-rho_uC_uL*a_unc1)))*binormal(a_unc1,a_L,rho_uC_uL)^(-1)		/*These expressions are the ones verified via Monte Carlo simulations*/
gen corr_C_L=normalden(a_L)   *(normal((1-rho_uC_uL^2)^(-0.5)*(a_unc1-rho_uC_uL*a_L)))*binormal(a_unc1,a_L,rho_uC_uL)^(-1)

***TABLE 10: RHO_{C,L}
reg u corr_C_0 corr_C_L if L==1 & uncensored_C==1,nocons
gen rcl = _b[corr_C_L]
 
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_L]						/*This is corr(u_C,u_L)*/
matrix semoments=semoments\_se[corr_C_L]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
											
												
biprobit (uncensored_C =$spec_U) (A =	$c_A),nolog	
predict a_unc2,xb1
predict a_App,xb2
scalar rho_uC_uA=e(rho)
scalar rho_unc_ap=e(rho)
	
gen corr_C_0_2=normalden(a_unc2)*(normal((1-rho_uC_uA^2)^(-0.5)*(a_App-rho_uC_uA*a_unc2)))*binormal(a_unc2,a_App,rho_uC_uA)^(-1)	/*These expressions are the ones verified via MC simulations*/
gen corr_C_App=normalden(a_App) *(normal((1-rho_uC_uA^2)^(-0.5)*(a_unc2-rho_uC_uA*a_App)))*binormal(a_unc2,a_App,rho_uC_uA)^(-1)

***TABLE 10: RHO_{C,A}
reg u corr_C_0_2 corr_C_App if A==1 & uncensored_C==1,nocons
gen rca = _b[corr_C_App]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_App]				/*This is corr(u_C,u_A)*/
matrix semoments=semoments\_se[corr_C_App]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

scalar rho_C_0_TBUB		=_b[corr_C_0_2]
scalar rho_C_App_TBUB	=_b[corr_C_App]

*gen Aw=1-R

biprobit (uncensored_C =$spec_U) (Aw =	$c_R),nolog	
predict a_unc3,xb1
predict a_Aw,xb2
scalar rho_unc_aw=e(rho)

egen c2=rmean(a_unc2 a_unc3)
gen  c3=A_App_3ml
gen  c4=A_Aw_3ml
scalar r23=rho_unc_ap
scalar r24=rho_unc_aw 
scalar r34=rho_aw_ap
scalar r34_2=(r34-r23*r24)/(sqrt(1-r23^2)*sqrt(1-r24^2))
scalar r24_3=(r24-r23*r34)/(sqrt(1-r23^2)*sqrt(1-r34^2))
scalar r23_4=(r23-r24*r34)/(sqrt(1-r24^2)*sqrt(1-r34^2))
matrix Sigma = (1, r23, r24 \ r23, 1, r34 \ r24, r34, 1)

mdraws, neq(3) dr(20) prefix(q) random
matrix C=cholesky(Sigma)
egen F3 = mvnp(c2 c3 c4), chol(C) dr(20) prefix(q) adoonly

gen corr_C_00=normalden(c2)*binormal((1-r23^2)^(-0.5)*(c3-r23*c2),(1-r24^2)^(-0.5)*(c4-r24*c2),r34_2)/F3
gen corr_C_ap=normalden(c3)*binormal((1-r23^2)^(-0.5)*(c2-r23*c3),(1-r34^2)^(-0.5)*(c4-r34*c3),r24_3)/F3
gen corr_C_aw=normalden(c4)*binormal((1-r24^2)^(-0.5)*(c2-r24*c4),(1-r34^2)^(-0.5)*(c3-r34*c4),r23_4)/F3

gen v=u-(rho_C_0_TBUB)*corr_C_00-(rho_C_App_TBUB)*corr_C_ap			/*use consistent estimates of corr.coeff. from other parts of the code*/

***TABLE 10: RHO_{C,R}
reg v corr_C_aw if Aw==1 & uncensored_C==1 & A==1,nocons
gen rcr = -_b[corr_C_aw]
scalar b_CA=_b[corr_C_aw]
scalar b_CAse=_se[corr_C_aw]

*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
matrix moments=moments\_b[corr_C_aw]
matrix semoments=semoments\_se[corr_C_aw]
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*
*&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*

matrix M=moments,semoments
matrix rownames M = 1_bF_V 2_bFv_V 3_corr_L_V 4_b_F_L 5_bF_A 6_corr_L_A 7_corr_R_A 8_corr_L_R 9_bF_C 10_sdu_C 11_EuC_L 12_EuC_A 13_EuC_Aw
matrix M=M[4..6,1..2]\M[1..3,1..2]\M[8..8,1..2]\M[9..11,1..2]\M[13..13,1..2]\M[12..12,1..2]\M[7..7,1..2]

preserve
	matrix list M
	svmat M
	keep M1 M2
	keep if M1!=.
	gen str20 namevar=""
	replace namevar="bF_L"            in 1
	replace namevar="bF_A"            in 2
	replace namevar="corr_L_A"        in 3
	replace namevar="bF_V" 			  in 4
	replace namevar="bFv_V"           in 5
	replace namevar="corr_L_V"        in 6
	replace namevar="corr_L_R"        in 7
	replace namevar="bF_C"            in 8
	replace namevar="sdu_C"           in 9
	replace namevar="EuC_L"           in 10
	replace namevar="EuC_Aw"          in 11
	replace namevar="EuC_A"   		  in 12
	replace namevar="corr_R_A"        in 13
	if order_expansion==2	{	
		save $data\moments_for_R_nonlin_expansion2,replace
	}
	else if order_expansion==3	{	
		save $data\moments_for_R_nonlin_expansion3,replace
	}
	
restore


***To get alpha
cap program drop get_alpha
program define get_alpha
	scalar sigma	=sqrt(1+s_csi^2)			/*This is [1+var(noise_SSA)] */
	gen xD=xbL_only*sqrt(1+s_omega^2+s_psi^2)+mu*female			
	ren a2 xA
	scalar rho_euA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)								/*We want the corr b/w epsilon and uA, not between uL and uA*/
	gen pDi=binormal(xA,xD,rho_euA)/normal(xA)
	gen m1i=normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,xD,rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,xD,rho_euA)^(-1)	
	gen m0i=-normalden(xD)*(normal((1-rho_euA^2)^(-0.5)*(xA-rho_euA*xD)))*binormal(xA,-xD,-rho_euA)^(-1) + ///
			rho_euA*normalden(xA)*(1-normal((1-rho_euA^2)^(-0.5)*(xD-rho_euA*xA)))*binormal(xA,-xD,-rho_euA)^(-1)	
	scalar sigma=sqrt(1+s_csi^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)
	su part1 if female==1 & A==1
	scalar p1_1=r(mean)
	su part1 if female==0 & A==1
	scalar p1_0=r(mean)
	scalar part1m=p1_1-p1_0
	su part2 if female==1 & A==1
	scalar p2_1=r(mean)
	su part2 if female==0 & A==1
	scalar p2_0=r(mean)
	scalar part2m=p2_1-p2_0
	su part3 if female==1 & A==1
	scalar part3m=r(mean)
	gen rho_FX		=part1+part2
	gen lambda_FX	=part3*female
	gen F=female
	probit R female college age racex2 experience married widowed doctor* walkra_R_1 dressa_R_1 stoopa_R_1 beda_R_1 bmi_R hosp_R i.year i.occ_interm bsx1-bsx9 agesq-occint8fem
	predict PI,pr
	gen psi=invnorm(PI)
	su psi 		if F==1 & R!=.
	scalar psi1=r(mean)
	su psi 		if F==0 & R!=.
	scalar psi0=r(mean)
	su rho 		if F==1 & R!=.
	scalar rho1=r(mean)
	su rho 		if F==0 & R!=.
	scalar rho0=r(mean)
	su lambda 	if F==1 & R!=.
	scalar lambda1=r(mean)
	scalar alpha1=((psi1-psi0)-(rho1-rho0))/lambda1		/*New alpha D-in-D: Feb 2024*/
	scalar list alpha1
end

cd $data

u V_pars,clear
sort draw
merge draw using R_pars
drop if _merge!=3
drop _merge
sort draw
merge draw using A_pars
drop if _merge!=3
drop _merge

egen nonconv=rsum(cvg_VL cvg_ALR cvg_C cvg_0L cvg_0A cvg_0Aw_R cvg_partcorr cvg_chol)
keep if nonconv==8

set seed 110968
gen uu=uniform()
sort uu
keep in 1/200

drop uu draw

ren alpha1 alpha_old
ren alpha1_alt alpha1		
su alpha1
scalar se_alpha1=r(sd)

corr bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L E_uC_R E_uC_A rho_A_R alpha1,cov 
matrix V_alpha=r(C)

corr bF_L bF_A rho_A_L bF_V bFV_V rho_V_L rho_L_R bF_C sigma_C E_uC_L E_uC_R E_uC_A rho_A_R,cov 
matrix V=r(C)
mat M3=vecdiag(V)'


***THIS IS EWMD (with the I matrix as weighting matrix)
if order_expansion==2	{	
	u $data\moments_for_R_nonlin_expansion2,clear
}
else if order_expansion==3	{	
	u $data\moments_for_R_nonlin_expansion3,clear
}
mkmat M1
mata: mata clear
mata:
void eval0(todo, b, v, g, H)
 {
		  x1 =(b[4]-b[5])/sqrt(1+b[1]^2+b[9]^2)										
		  x2 =(b[4]-b[5]-b[6])/sqrt(1+b[1]^2+b[9]^2+b[2]^2)							
		  x3 =(1+b[1]^2+b[9]^2)/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2+b[2]^2))	
		  x4 =-(b[5]-b[7])/sqrt(1+b[1]^2+b[9]^2)											
		  x5 =(b[4])/sqrt(1+b[1]^2+b[9]^2)											
		  x6 =b[9]^2/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[1]^2+b[9]^2))					
		  x7 =-1/(sqrt(1+b[1]^2+b[9]^2)*sqrt(1+b[11]^2))								
		  x8 =b[8]*b[4]+b[10]															
		  x9 =sqrt(b[8]^2+b[3]^2)														
		  x10=b[8]/sqrt(1+b[1]^2+b[9]^2)												
		  x11=b[8]/sqrt(1+b[11]^2)
		  x12=b[8]/sqrt(1+b[1]^2+b[9]^2+b[2]^2)
		  x13=-1/(sqrt(1+b[1]^2+b[9]^2+b[2]^2)*sqrt(1+b[11]^2))
		  M1=st_matrix("M1")
		  M2=st_matrix("M3")
		  m1 =M1[1,1]
		  m2 =M1[2,1] 
		  m3 =M1[3,1]
		  m4 =M1[4,1]
		  m5 =M1[5,1] 
		  m6 =M1[6,1]
		  m7 =M1[7,1] 
		  m8 =M1[8,1]
		  m9 =M1[9,1]
		  m10=M1[10,1]
		  m11=M1[11,1]
		  m12=M1[12,1]
		  m13=M1[13,1]
		  v1 =M2[1,1]
		  v2 =M2[2,1] 
		  v3 =M2[3,1]
		  v4 =M2[4,1]
		  v5 =M2[5,1]
		  v6 =M2[6,1]
		  v7 =M2[7,1] 
		  v8 =M2[8,1]
		  v9 =M2[9,1]
		  v10=M2[10,1]
		  v11=M2[11,1]
		  v12=M2[12,1]
		  v13=M2[13,1]
	v=((x1-m1)^2/1)+((x2-m2)^2/1)+((x3-m3)^2/1)+((x4-m4)^2/1)+((x5-m5)^2/1)+((x6-m6)^2/1)+((x7-m7)^2/1)+((x8-m8)^2/1)+((x9-m9)^2/1)+((x10-m10)^2/1)+((x11-m11)^2/1)+((x12-m12)^2/1)+((x13-m13)^2/1)
 }
S = optimize_init()
optimize_init_evaluator(S, &eval0())
optimize_init_evaluatortype(S, "v0")
optimize_init_params(S, (1,1,1,1,1,1,1,1,1,1,1))
optimize_init_which(S, "min")
optimize_init_which(S, "min")
optimize_init_conv_ptol(S, 1e-16)
optimize_init_conv_vtol(S, 1e-16)
b = optimize(S)
b'
st_matrix("beta_ewmd",b')
end


scalar mu		=beta_ewmd[4,1]
scalar s_omega	=beta_ewmd[1,1]
scalar s_psi	=beta_ewmd[9,1]	
scalar s_v		=beta_ewmd[2,1]
scalar s_csi	=beta_ewmd[11,1]	
u data_to_est_alpha_nonlin,clear
get_alpha
matrix beta=beta_ewmd\alpha1
svmat beta
keep beta
drop if beta==.
ren beta1 beta
gen str20 namepar=""
replace namepar="sigma_omega"      in 1
replace namepar="sigma_v"   	   in 2
replace namepar="sigma_phi"        in 3
replace namepar="pi"  		       in 4
replace namepar="gamma"  	       in 5
replace namepar="tau"   	       in 6
replace namepar="theta"  	       in 7
replace namepar="lamda"  	       in 8
replace namepar="sigma_psi"        in 9
replace namepar="psi"   	       in 10
replace namepar="sigma_csi"        in 11
replace namepar="alpha"        	   in 12
order namepar beta

*****************************
if order_expansion==2	{	
	save results_nonlin_expansion2,replace
}
else if order_expansion==3	{	
	save results_nonlin_expansion3,replace
}


**************************************************************
***This is Table OA6******************************************
**************************************************************
u $data\moments_for_R,replace
sort namevar
ren M1 M1e1
drop M2
save $data\mt1,replace
u $data\moments_for_R_nonlin_expansion2
sort namevar
ren M1 M1e2
drop M2
save $data\mt2,replace
u $data\moments_for_R_nonlin_expansion3
sort namevar
ren M1 M1e3
drop M2
merge namevar using mt1
drop _merge
sort namevar
merge namevar using mt2
drop _merge
order namevar M1e1 M1e2 M1e3
erase mt1.dta
erase mt2.dta
mkmat M1e1 M1e2 M1e3,mat(mom) rownames(namevar) 

esttab matrix(mom, fmt(%9.3f)) using "$figures/tableOA6_1.tex", ///
					title("Moments") collabels("Baseline" "Quadratic" "Cubic") replace ///
					postfoot(\hline\hline \end{tabular} \end{table} )

u $data\struct_pars_estimates_ewmd,replace
sort namepar
ren beta be1
save $data\mt1,replace
u $data\results_nonlin_expansion2
sort namepar
ren beta be2
save $data\mt2,replace
u $data\results_nonlin_expansion3
sort namepar
ren beta be3
merge namepar using mt1
drop _merge
sort namepar
merge namepar using mt2
drop _merge
order namepar be1 be2 be3
erase mt1.dta
erase mt2.dta
mkmat be1 be2 be3,mat(est) rownames(namepar) 

esttab matrix(est, fmt(%9.3f)) using "$figures/tableOA6_2.tex", ///
					title("Parameters") collabels("Baseline" "Quadratic" "Cubic") replace ///
					postfoot(\hline\hline \end{tabular} \end{table} )

